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^ Abstract 

on' 

We study general stochastic birth and death processes including delay. We develop several ap- 
proaches for the analytical treatment of these non-Markovian systems, valid, not only for constant 
delays, but also for stochastic delays with arbitrary probability distributions. The interplay be- 
tween stochasticity and delay and, in particular, the effects of delay in the fluctuations and time 
correlations are discussed. 
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I. INTRODUCTION 



Stochastic modeling plays an important role in many areas of science, such as physics, 
ecology or chemistry Stochasticity may appear due to the lack of complete knowledge 
about all the relevant variables, the precise dynamics of the system or the interactions with 
the environment. In some cases, one can obtain a compact description of a complicated 
system considering only a few relevant variables but at the expense of losing deterministic 
predictability. Often, probabilities for some fundamental processes can be assigned on the 
basis of symmetries and other considerations, or on empirical analyis, and the dynamics of 
the process can be derived bottom-up. 

Stochasticity appears together with delay terms in many situations of interest, such as 
gene regulation (2 - 4 ] , physiological processes jlj] or postural control 

a a. 

The combined 

effects of stochasticity and delay are, however, not completely understood. From the math- 
ematical point of view, stochastic processes including delay are difficult to analyze due to 
the non-Markovian character. Most of the previous approaches have focused on stochastic 



differential equations, that consider continuous variables [8Nl2|. or random walks in discrete 
time [13I [ijj], where delay can be taken into account increasing the number of variables. 
Models with discrete variables but continuous time are the natural description of many sys- 
tems such as chemical reactions, population dynamics or epidemic spreading. In some cases, 
discreteness can be a mayor source of fluctuations, not well captured by continuous models 
15| . The approach with discrete variables and continuous time was used in 4|, II6MI8I] . Most 
often, the delay time is taken to be a constant with zero fluctuations. This is not very 
realistic in the applications, since it is unusual to have a deterministic delay when the rest 
of the dynamics is stochastic. We will take this consideration into account by allowing the 
delay times to be random variables with arbitrary probability density functions. 

In this work we study some simple, yet general, stochastic birth and death processes 
including delay. We will develop tree different approaches to the analytical study of this kind 
of non-Markovian processes, in the general case of stochastically distributed delay: a direct 
approach in subsection (III AD . an effective Markovian reduction in subsections (III Bp and 
(III CI) , and a master equation approach, together with a time- reversal invariance assumption, 
in section (IIIII) . The first direct approach method is interesting for its simplicity, but its 
application is limited to systems with first order reactions and without feedback. The second 
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one, effective Markovian reduction, is rather ffexible and general and its development is one 
of the main advances of this paper. The last master equation approach complements the 
previous, giving information about the full probability distribution. The main limitation 
of all the approaches is the need to assume that completion times for delayed reactions 
are independent random variables (independent of each other and of other variables of the 
system), although the initiation rates may depend on the state of the system, allowing, 
for example, for feedback and crowding effects, so we do not consider this limitation to be 
very relevant for practical applications. Although our methodology is rather general, we 
present it here using specific examples that have been grouped in two categories: delay in 
the degradation (section [XT]) and delay in the creation (section HTT1) . We end the paper with 
a brief discussion and comments in section IIVI Some more technical details are left for the 
two appendices. 

II. DELAYED DEGRADATION 

We will start by studying simple stochastic birth and death processes that include delay 
in the degradation step. A process of this type was proposed in [4] as a model for protein 
level dynamics with a complex degradation pathway. 

A. Simple Case 

We consider first the simplest possible process including delayed degradation: 

c 

0^X, X^0, (1) 
r 

that is, a particle X is created at a rate C and disappears ("dies" or "degrades") a time 
t after created. We allow the delay time r to be randomly distributed i.e. the lifetimes r 
of the created particles are random variables, that for simplicity we consider independent 
and identically distributed, with probability density fij). Although not considered in this 
paper, the case of non-identically distributed delay times, in particular a probability density 
that depends on the time from birth, can also be treated. However, as commented above, 
the case of non-independent delay times does not seem to be tractable with the methods we 
present below. 



3 



We note first that distributed delay is completely equivalent to degradation at a rate that 
depends on the "age" a (time form creation) of the particle, i.e., processes 

7(a) 

X^0, and X^0, (2) 

T 

are equivalent if the rate 7(a) and the probability density of the delay f(r) are related by: 



7(aJ 



F(a) 



f(r)= 1 (r)e-fo da ^ a \ (3) 



with F(t) = 1 — F{t) being F(t) = Prob(r < t) = J" r drf{r) the cumulative distribution 
of the delay-time. This is so because 7(a)da is the probability of dying at the time interval 
(a, a + da), if the particle is still present at a, and so it is nothing but the probability f(a)da 
that the delay time r belongs to that same interval conditioned to the particle still being 



alive at time a, an event with probability F(a). In the notation of 19j, 7(a) is nothing 
but the conditional failure rate. We take t = as the time origin, so the number of alive 
particles at time t is n(t) = for t < 0. Let P(n, t) the probability of n particles being alive 
at time t. In the remaining of this subsection we assume that there is no feedback, in the 
sense that the creation rate C is independent on the number of particles n, but, for the sake 
of generality, we do allow it to be a function of time C(t). The non-feedback assumption 
allows us to obtain a full analytical solution. As shown in the appendix, independently of 
the form of the delay distribution, P(n, t) follows a Poisson distribution 

P(n,t) = e-W))W$£ l (4) 

with average (n(t)) = f* dt'C(t')F(t — t'). If the creation rate, C(t), is independent of time, 
a steady state is reached, in which the average number of particles is (n) s t = C(r), again 
independently of the form of the delay distribution. 

We will now compute the time correlation function. We shall see that its analytical 
expression does depend on the form of the delay distribution. We start from the relation: 

(n(t + T)\n(t)) = (n new (t + T)\n(t)) + (n old (t + T)\n{t)), (5) 

with n new {n id) particles created after (before) t. n new can be computed exactly as before 
(now taking t as the time origin), so we have: 

(n new {t + T)\n(t))= [ dt'C(t + t')F(T-t'). (6) 
Jo 



The evolution of the number of particles already present at t depends on the age a of these 
particles. Their survival probability until time t + T can be written as: 

P(alive at t + T|alive at t) = / daP(&ge = a|alive at t)P (lifetime > a + T|lifetime > a) 

Jo 

' da ~ + T "> - & *eCtf)F(t +T-?) 



o 



f* dt'C(t')F(t - t>) F{a) C dt'C{t')F{t - V) 



where we used P(a\b) = P p$ , so we find: 

/ *a.ru m\ t Jldt'C{t')F{t + T-t') 

{Tioidit + T)\n(t)) = n[t)— — 7 . 8 

V 1 M K " KJ f*dt'C(t')F(t-t>) 1 ; 

From this, one easily obtains the correlation function: K[n](t,T) = {n(t){n(t + T)\n(t))) — 
(n(t))(n(t + T)) 

K[n){t,T) = [ dt'C(t')F(t + T-t'), (9) 
Jo 

If C(t) = C, independent of time, a steady-state can be reached with correlation function 
K st [n](T) = lim^oo K[n](t, T). For a constant rate 7, which would be equivalent to an 
exponential delay distribution /(r) = 7e _7r , it has the usual exponential decay K st [n](T) = 
(C/7)e~ 7T . For a fixed delay time r , corresponding to /(r) = 8(r — r ), the correlation 
function is a straight line K st [n](T) = C(r — T) for T < r and K st [n](T) = for T > r . 
For other distributions of delay time, the correlation function adopts different forms, but it 
is always monotonically decreasing. In figure (CQ) we plot the correlation function for two 
different types of distribution of delay, for different values of the variance of the delay. We 
see that the distribution with fatter tail displays a slower asymptotic decay, and that the 
decay is slower as the variance of the delay increases. Numerical simulations, performed with 
a conveniently modified version of the Gillespie algorithm [20J], are in perfect agreement with 
this exact result, providing a check of its correctness. We remark that the functional form 
of the decay of the correlation function depends on the delay distributed and can differ from 
the exponential decay found in systems without delay. 



B. More elaborated Case 

We now consider a process including both instantaneous and delayed degradation steps: 

c 7 D 

0^X, X^0, =>0, (10) 




t 



FIG. 1: Steady state correlation function, Eq.Q, as a function of time, plotted in logarithmic 
scale, for two different types of delay distribution, gamma and lognormal, for two values of the 
variance of the delay: = 0.2 (left panel) and a\ = 5 (right panel); in both cases the average 
delay is (r) = 1 and the creation rate is C = 1. We also plot a exponential decay with exponent 
one (dot-dashed line), for comparison. Note that delay distributions with larger variance and fatter 
tayls display slower asymptotic decay. (Online version in colour.) 

this is, particles are created at a rate C and each particle can be eliminated by two processes: 
i) instantaneous degradation at a rate 7; ii) delayed degradation, initiated at a rate D but 
completed only a time r after initiation. Again, we will allow the delay-degradations times 
to be random variables that, for simplicity, will be independent and identically distributed 
with probability density function /(r). 

For the process to be completely defined, one has to specify if a particle that initiates 
delayed-degradation at time t and thus will disappear at t + r (this kind of particles will 
be called "infected"), can also disappear before the completion of this reaction, through 
instantaneous degradation. In the most general case, this can happen at a rate 7', not 
necessarily equal to 7. Note that, in the case of first-order degradation (7' not dependent 
on the number of particles n), this instantaneous degradation is completely equivalent to a 
system with 7' = 0, after modifying the distribution of the delayed-degradation times in the 
following way: 

/( T ) _> e -^/(r) + e-^iF{r). (11) 

That is, when instantaneous degradation is added to infected particles, the probability that 
the lifetime is equal to r has two contributions: (i) a particle initially has a lifetime r 
(probability density /(r)) and survives up to this time (an event with probability e~ 7 ' r ); (ii) 
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a particle has a lifetime larger than r (probability F(t)), but survives up to r (probability 
e _7 ' r ) and then undergoes instantaneous degradation (at rate 7'). The consideration of 
these two contributions leads straightforwardly to Eq.f lllj) . We see that omitting first order 
instantaneous degradation of infected particles comprises no loss of generality, given that 
the treatment is valid for general distributions of delay. 

If D and 7 are independent of n, the process is equivalent to the one-variable system 
discussed in the previous subsection III Ap with a conveniently modified distribution of delay: 

f(r) -> e - (7+D)r 7 + f dt'e^ +D)t 'Df{T - t'). (12) 
Jo 

This comes from the fact that a particle may disapear at time r because it did not disapeared 
or was infected before and is degraded instantaneously (probability density e _< - 7+D - )T 7) or 
because it got infected at some previous time (£') with an appropriate lifetime (r — t', 
probability density JJ" dt'e~(~ /+D ^ t ' Df(r — t')). This includes as particular cases the ones 
studied in [17|, [22|. The results of subsection ( III Al) allows us to obtain the full solution also 
in the general case of distributed delay. If D or 7 depend on n the processes are not anymore 
equivalent, two variables are necessary and a new approach is needed for the analysis. In 
the following we develop this method. We will also consider the case in which the creation 
rate C depends on the number of particles. 

The full process corresponds to the following two-variable system: 

c 7 D 

— X A , X A — 0, X A — Xj + Z, Xj => 0, (13) 

T 

where we have split the proteins into two types: Xi are infected particles that will die 
precisely at a time r (itself a stochastic variable) after being infected and X A are non- 
infected ("active") particles (so X = X A Ulj). We allow the rates to depend on n A , the 
number of X A , active, particles, but not on nj, the number of Xj, infected, particles which 
are considered to be "inert" ; this condition will be realxed in the next subsetion. Following 



17j |. we have introduced the auxiliary particles Z whose number is given by the stochastic 
variable nz(t). The introduction of Z will allow us to obtain the properties of rii by using 
the relation: 

nj(t) = £ dt ,d ^ls{t',t), (14) 

where the discrete process nzit) is a sequence of step (Heaviside) functions and its derivative 
must be understood as a series of Dirac-delta functions. Here we have introduced the family 



of "survival" stochastic processes s(t',t) defined in the following way: first, for each t' 
we obtain a value of r(t') independently drawn from the distribution /(r). Next, we set 
s(t',t) — 1, if t G (i', t' + r(t / )) jL and s(t',t) = 0, otherwise. This can be considered as the 
indicator function of a virtual [24J particle that is infected at t' and survives up to a time 
t' + r(t'). It follows from the definition that: 

(s(t ls t)) = F(t-ti), (15) 

(s(M)s(* 2 ,0) = <( (16) 

(s(ti, max{t, i'})) if t\ = t 2 , 



Expressions ( JT41T16]) are the main advances of this section and provide us with the neces- 
sary tools to derive the main properties of the stochastic process ( ITU . In the case considered 



m 



171 ] there is a fixed delay (/(r) = 5{r —To)) and no instantaneous degradation of infected 
particles (7' = 0), so one has simply nj(t) = nz{t) — nz{t — r). The inclusion of the survival 
process s(t', t) allows us to consider the general case of distributed delay and rates depending 
on the state of the system. 

Note that the process followed by {n A , nz} is Markovian as the delay only appears in 
variable Xj, so the properties of nz can be obtained using Markovian methods, and the 
properties of the variable nj can be derived afterwards using ( Ti4iri6l) . In particular, the first 
moments follow: 



(n/(t)> = l t Jt' d ^^{s(t'M (17) 



(n/(t)n/(t + r)) = / dh I dt 2 d {nZ ^ t f h)) (s(t u t)s(t 2 ,t + T).) (18) 



—00 j —00 



Using standard Markovian methods [l], one can prove that the process {n^, %} is described 
by the master equation: 

dP{nA dt nZ,t) = (E A l - l)C(n A )P(n A ,n z ,t) + (E A -l)-y(n A )P(n A ,n z ,t) 

+ (E A E z l )D(n A )P(n A ,n z ,t), (19) 

with Ei the step operator, £^/(n», %) = /(n* + 1, nj). In this section, we allow the creation 
rate C to depend on the number of X^-particles, constituting a feedback term on the number 
of "active" particles. From the master equation one easily derives the equations for the 
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moments, the first of them read: 
d(n A ) 



dt 
d{n z ) 

dt 
d(nl) 

dt 
d(n%) 

dt 

d{n A n z ) 
dt 



(C(n A ))-((j + D)n A ), 
(Dn A ), 

2(2(n A + l)C(n A )) - ((2n A - l)n A ( 7 + D))), 
2(Dn A n z ) + (Dn A ) 

{C(n A )n z } - ((7 + D)n A n z ) + (D(n 2 A - n A ))). 



(20) 
(21) 
(22) 
(23) 
(24) 



In the case that C(n A ) is a linear function of n A and 7 and D do not depend on n A (and 
none of them depend on rii or n z ), the system of equations is closed and can be solved. For 
non-linear systems, we will make use of van Kampen's expansion This is a standard 
systematic expansion of the master equation, that consists on assuming a deterministic, and 
a stochastic part for the variables that scale differently with a large parameter, Q (typically 
the volume or system size) i.e. n A = Q<p A (t) + Q 1 ^ 2 ^, n z = fl(p z (t) + VL 1 / 2 ^ Z . One can then 
write the master equation for the new variables £,4, ^ z and expand in powers of Vt~ l l 2 . The 
method is generically valid, provided that the rates depend on the variables only trough 
n a /fl (plus higher orders in fi -1 ; a common factor depending on Q multiplying all rates is 
also acceptable), which is fulfilled by most systems of interest, and that the macroscopic 
equations have a steady state as a single attractor. The equations for the macroscopic 
components are: 



dt 
d(f> z 
dt 



^(0A)-[7(0A)+^(0A)]0A, 

DU A U A . 



The stochastic contributions, to first order in f2 1//2 , read: 
d(U) 



dt 
d(Zz) 

dt 
d(ti) 

dt 

d(? z ) 
dt 
d(UZz) 
dt 



- 7 + J D-C"(0 A )j (U), 

-2 [7 + D - C'{(j> A )] (e A ) + (i + D)<p A + C{4> A ), 
2D{UZz)+D4> A , 

7 + D - C'(<f> A )] + D((e A ) - 0a), 

9 



(25) 
(26) 

(27) 
(28) 
(29) 
(30) 
(31) 



with D = D{4>a) + D'((f)A)(pA, 7 = 7(^.4) + i'( ( Pa)4 > a- Usually, for the ansatz about the 
scaling of the variables to work (and so the expansion), the equations for the macroscopic 
components must have a single stable fixed point. In this case, however, the equation for 
2 does not have a fixed point, and (f> z (t) and (£§(i)) g row without bound. This grow, 
nevertheless, is consistent with Js^^lt) = ^(^°)' anc ^ the expansion can still be applied. 

( T2~7ll3~Tj) is a system of closed linear equations and so can always be solved. To compute 
the time correlations of nj from Eq. (1181) we need the time correlations of n z . We note that: 

(n z (ti)n z (t 2 )) = nzinz2P(nz2,t 2 ;n Z i,t 1 ) (32) 

= n zin Z2 P{n Z2 ,t 2 \n Zl ,n A ,t 1 )P{n ZA ,n A ,t 1 ) 
= ((n z (t 2 )|nz(ti),TiA(ti)>n z (t 1 )> , (33) 



and that (nz(t 2 )\nz(ti), n A (ti)) (for t 2 > £1) can be obtained integrating (I20H21I) or 
In the general, non-linear, case, using first order van Kampen's expansion, one obtains, over 
the steady state: 



{nz(ti)nz(t 2 )) = ^ 2 Mti)4>z(h) + « 



u 



, (34) 



with u = 7 + D — C'{4>A,st) an d 4>A,st the solution of C(4> A ) — (7 + P j )4 ) a- The derivative 
that appears in ( ITS]) is: 

rf 2 (n z (ti)n z (t 2 ))st 



« 2 



^n^z)^ - " 1 * 1- * 8 ' + D^ AtSt 6(ti - t 2 ) , (35) 



with {^Aiz)st = D4>A,st— — ^2v? ^ A ' 8t ■ Putting all the pieces together, one finally obtains: 



^ s tN(t) = (n / (t )n/(to + t)) s t-(n/)f t (36) 

roc POO POO 

= nD<p AM / dt'F(t + t') + nDu(^ z ) st / / drF(s)F(r)e- u \ t+s - r \. 
Jo Jo Jo 

Proceeding in a similar way, one can derive: 

POO 

K st [n A ,nj}(t) = niU&BtU / dt'e- u ^F(t') (37) 

Jo 

poo rt 

K st [nj,n A ](t) = Sl(UZz)«u / dt' e~ ut ' F{t + t') + SlD^st / df 'e~ ut ' ' F(t - t'), (38) 

Jo Jo 

K st [n A )(t) = n(Q at e- ut , (39) 
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with K st [n u , n v ](t) = (n u (t + t)n v (to)) st — (n u ) st (n v ) st . This finally allows to express the 
correlation function for the total number of particles, n = n A + ni, as: 

K Bt [n](t) = K st [nj](t) + K st [n A ,m](t) + K st [n h n A }{t) + K st [n A }(t). (40) 

In this case, the average of n again depends only on the average delay, (n) st = ^70^(1 + 
D(r)), but the second moment depends on the delay distribution in a more complicated 
way, through factors involving the integral of F(t). 

In figure ([2]) this result is compared with numerical simulations, showing a very good 
agreement. Note that the treatment of the delayed reactions is exact, the only approximation 
coming from the use of van Kampen's expansion, which is needed when non-linearities are 
present, but whose error scales as Vt~ 1 / 2 . Like in the previous case, the process in which the 
distribution of delay has fatter tail shows slower decay for the correlation function. 



g 2 = 0.2 a 2 = 5 

X T 




2 4 5 10 15 20 



t t 

FIG. 2: Steady state correlation function for the total number of particles as a function of time, 
plotted in logarithmic scale, for two different types of delay distribution, gamma and lognormal, 
for two values of the variance of the delay: 0.2 (left panel) and 5 (right panel); in both cases 
the average delay is (t) = 1. The insets show the time correlation for the number of "infected" 
particles, Xj, which gives the largest contribution to the difference between different distributions. 
Symbols come from numerical simulations and lines from the theoretical analysis Eqs. (|37ti^0|) . The 
creation rate is C(n A ) = 1+ ^ A ^ ' pa ram eters values are: O = 100, c = 1, e = 0.4 and D = 7 = 1. 
(Online version in colour.) 
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C. Full feedback 



X, X^0, (41) 



We now consider the case in which the creation rate depends on all present particles 

C(n) 

X, X 

T 

with n the total (inert+active) number of X-particles. As noted before, this single-variable 
model can account for instantaneous plus delayed degradation, in the case that the degra- 
dation and "contagion" rates, 7 and D before, do not depend on the state of the system. 
For simplicity, we restrict our attention to this case. This process can be treated with the 
approach of the previous subsection introducing the additional variable Z, 

C{n) 

®^X + Z, X^0, (42) 

T 

with riz(t) the corresponding random variable giving the number of Z particles. We see 
that: 

n(t) = f dt'^jPs(t',t), (43) 

with s(t',t) the same as in the previous section. The probability distribution for {n,riz} 
follows a master equation of the form: 

dP ^ n ^ = (E-'Ez 1 - l)C(n)P(n z , n, t) + (E - l)g(n z , n)P(n z , n, t). (44) 

Details of the derivation of the master equation in systems with delay are given in the 
appendix. Here, g(n z ,n) = J °° dt' (C(n(t — t'))\n z (t),n(t)) f(t'), with f(t) the probability 
density of the delay distribution, although, since we are only interested in the properties of 
variable n, we will not be using this expression. The key step in this case is to note that 
Eq. fj44l) allows us to derive the statistical properties (moments and correlations) of n z (t) as 
a function of those of n(t). Then, using (143]) we will be able to self-consistently derive the 
properties of n. More specifically, the approach proceeds as follows: 

Summing Eq. (l44j) over n, we can obtain an equation for the evolution of P(n z ,t), but that 
still depends on n (in this step the contribution of the second term in Eq. (144|) vanishes): 

dP{ ^ ,t] = ( E z l ~ llE^K.M) = (E, 1 - l){C(n(t))\n z ,t)P(n z ,t). (45) 

n 

The two times probability distribution P(n Zl ,ti;nz 2 ,t 2 ) follows a similar equation. Condi- 
tioning carefully, summing over the variable n and considering separately the case t\ = t 2 

12 



(which turns out to be singular), we find: 

d2P{nZ dtldt7 2M) = C^ 1 - i)^ 1 - l)<C(^(ti))C(^(t 2 ))|^ 1 ,t 1 ,^ 2 ,t 2 > x 

P(n Zl ,t 1 ;n Z2 ,t 2 )+ (46) 
5{h - t 2 ) [{I - E Z2 )5 nZi , nz ,E z l + (1 - E z l)5 nz ^ nz2 ] (C{n A )\n Zl )P{n Zl M) 

From f l4"3]) we easily obtain: 

(n(t)> = j^^^is^t)) (47) 

(n(t)n(t')) = f dh f dt / 2 M^fe)) (s(tl|t)s(t2|t0) , (48) 



(C(n(t))), (49) 
(C(n(* 1 ))C7(7i(* 2 ))) + 5(t! - t 2 ) (C(r^i))) (50) 



While (pi |46D imply: 

eft 

^(nz(ti)nz(^)) 
(iii(it2 

And we finally obtain the following set of integral equations for the moments: 

(n(t)) = f dt^Cinit^Fit-h) (51) 

J — oo 

t rt' 



(n(t)n(t')) = / dt x \ dt^C{n{t x ))C{n{t 2 )))F(t-t x )F(t! - 1 2 ) 

J —oo J — oo 

+ f dti(C(n(ti)))F(max{t,t'}-ti), (52) 

J — oo 

In the case of linear feedback, C(n) = a + bn, this system of equations is closed. For non- 
linear systems, one can use van Kampen's expansion as explained above. In the steady state, 
one finds: 

(n) st = O0 st , <f>st = C(0 st )(r) (53) 

/•oo r ft+x poo 

K st [n](t) = dx dyK st [n}(t + x-y)F(x)F(y) + / ^ st [n](t + z - y)F(x)F(y) 

JO Uo Jt+x 

POO 

+ fiC(&t) / dif(Hi) (54) 



o 



Eq. ( 153]) shows that the steady state number of particles depends only on the average 
delay. Eq. ( 154"]) shows that the correlations depend on the delay distribution in a non-trivial 
way. The analysis of this equation is left for future work. 
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III. DELAYED CREATION 



We now turn our attention to the case in which the creation reaction, that is initiated 
stochastically, takes a finite time to be completed. For simplicity we assume that the degra- 
dation reaction is instantaneous. Schematically, we have: 

C(n) 7 



X, X 



T 



(55) 



In this case, if the creation rate does not depend on the number of particles, n, then the delay 
in the creation is completely irrelevant, since the probability that a new particle appears at 
time t is equal to the probability that its creation started at a time t — r, but this equal to 
the probability that a particle starts its creation at time t (with a shift in the time if C is 



time-dependent 



so the process is completely equal to one with instantaneous creation. 



Following 18[ we will adopt here an approach different from that of the previous sections, 



that, besides the moments, will allow us to obtain an expression for the full probability 



distribution. For completeness we will exp 
considerations, the reader is referred to 
equation of the process (I55p is: 



ain here the method in some detail. For additional 



18]. In the appenix it is shown that the master 



(E-l)[ 7 nP(n,t)] + (E- 1 -l) 



dP{n, t) 
dt 

The master equation ( 1561) can be written as: 
dP(n,t) 



V / drC(n')P(n',t-T;n,t)f(T) 
n'=0 Jo 



(56) 



dt 



(E-l) [ 7 nP(n, t)} + (E- 1 - 1) [C(n, t)P(n, t)] 



(57) 



where the effective creation rate, C(n,t), is given by: 



C(n,t) 



drf(r)(C(n'(t-r))\n(t)). 



(58) 



The conditional probability P(n, t\no, to) follows a master equation identical to ( l56l) with 
all the probabilities conditioned to no at time to. From it, and using that (n(t)\n(to)) = 
nP(n, t\n(to),to), we obtain the following evolution equation for the conditional average: 

d(n(t)\n(t )) 



dt 



7 <n(f)|n(*o)>+ / drf(r)(C(n(t-T))\n(t )), 



(59) 



for t > 0, with initial condition (n(tQ)\n(t )) = n(t ). 
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The knowledge of the steady value C 8 t(n) = lim t ^ 00 (C(n'(t — r))\n, t)) = (C(n f ), — r\n) st , 
allows the calculation of the steady-state probabilities P st {n), obtained by imposing d p 0M) — 
in Eq.(l5T|h as 



hi 



pm = a.(o) n = ^ n e« (*>, (60) 

fe=0 ' V 7 ' k=0 

P st (0) is fixed by the normalization condition. All is left to do now is to compute the effective 
creation rate C s t(n). 

The effective creation rate will be computed using expression fl59l) . In the general case of 
nonlinear creation rate, we will use van Kampen's expansion to linearize C(n) around the 
macroscopic component of n. We have: C(n) = QC((f)) + f2 1,/2 C"(0)£, so 

(C(n'(t - r))\n(t)) = QC{<f>{t - r)) + ^C'm - r))(£'(t - r)|£(f)> (61) 



/■oo 

- 7 0(t)+ / drf(r)C^(t-r)), (62) 

POO 

-liat'Mt)) + / drf(r)C m - r)) (&) ~ r|^(t)> (63) 



using (J59|) we obtain: 

rf0(t) 

dt 

dt' 



o 



Equation (|62|) is in general a non-linear integro-differential equation, that can be difficult 
to solve. Here, however, we will focus on the cases in which fl62l) has a stable steady state 
as a single attractor, which is the solution of 70 = C{4>). This is the regimen in which the 
validity of van Kampen's expansion is guaranteed. 

We reach now a delicate point. Eq. (|63l is a (linear) integro-differential equation. To 
solve it, we would need an initial condition in the whole interval (— oo,t) but we only know 
a one-time condition (£(£' = = £(£)• We will circumvent this difficulty by assuming 

that, over the steady state, the system is statistically invariant under time- inversion, which 
implies (£(t + = — ^i) !£(£))• This condition, together with the value of £ at 

time t, allows to find the solution of ( 1631) . The time-reversal invariance assumption in the 
steady state is fulfilled by any Markovian system that follows detailed balance. Our system 
follows detailed balance (as any one-step process |l|), but, due to the presence of delay, it 
is not Markovian. So the time-reversal invariance is an assumption, whose validity needs to 
be checked. It was shown in [3] that in this system the assumption is approximately valid. 
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In the case of constant delay, f(t) = 5(t — r), the time-reversal symmetric solution of 



§B is mm- 



(£(t + ti)|£(t)> = Z(t)h(t 



r 



h{h) = < 



e-^- k ^h(kr) - a f£ dt' h(t' - r)^'-^, if kr < t x < (k + l)r, k = 1, 



A = ^ff^\ ( = ' L ~^, " = -C'{<j> s t) 

a 

and using (|6T|) we finally obtain: 

C(n) = [C(0 st ) - st C / (0 st )/i(r)] + C'(<f> at )h(T)n, (66) 

where si is the steady state solution of fl62|) . From Eq.f l60|) one can obtain the steady-state 
probabilities P st (n). The mean value and variance are given by: 



(n) st = Q(j) s t (67) 

(68) 



„2 ( n )l 



From fl68|) one can see that, interestingly, in the case of negative feedback (C"(0 st ) < 0), 
as the delay is increased the fluctuations change from sub-Poissonian (a 2 < (n)) to super- 
Poissonian (cr 2 > (n)). This is illustrated in figure (EJ), where we also show the line of the 
Hopf bifurcation for the deterministic system. It is usually obtained that a negative feedback 
reduces the magnitude of the fluctuations [23], when delay present we see that this negative 
feedback can change totally its effect, giving rise to an increase of the fluctuations. 
The time correlation function can also be obtained from ( 164]) . as: 

K st [n}(t) = a 2 st h(t) (69) 

In the case of negative feedback it becomes non-monotonic, developing peaks of alternating 
sign at approximately multiples of the delay, signaling the presence of stochastic oscillations. 
For positive feedback, the time correlation is always positive, but not necessarily monotonic. 

We will finish by noting that the "effective Markovian reduction" method used in the 
previous section can also be used for the case of delay in the creation with feedback. To be 
completely general, we allow two delays, one in the creation (with probability density f c {t)), 
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and one in the degradation (with probability density fd(t))- The process is schematized as 
follows: 

C(n) 

-> =*X, X^0, (70) 

with r c / d random variables distributed according to f c /d(t). With the addition of two new 
variables, the process can be rewritten as: 

C(n) 

(&^Z + Y, F^X, X^0, (71) 
t c r d 



which allows us to note that: 



n(t)= f dt'^±~s(t',t). (72) 



— oo 



In this case, the survival function s(t',t) is defined as: s(t',t) — 1, if t G (£' + T c (t'),t' + 
T c(t') + Td(t')) , and s(t',t) = 0, otherwise, being t c (£') and 7d(t') random times obtained 
from the corresponding pdf's / c (r c ) and fd(T d ). s(t',t) is equal to one if a virtual particle 
that initiated its creation at time t' finished it at some intermediate time t" < t and since 
then had a lifetime greater that t — t", so that it is still alive at t, being zero otherwise. It 
follows that: 

(3(*i,*)> = f ^ dt'f^F^t-h-t') (73) 
Jo 

(s(t u t)s(t 2 ,t')) = { (74) 
[ f^ om } - h dt"f c (t")F d (msx{t,t'} -t x - t'% if ^ = t 2 . 

In the case that the creation rate C(n) does not depend on the number of X-particles, the 
number of Z-particles follows a Markovian process (Poisson process), and the properties of 
n can be derived from (1721) . If the creation rate depends on the number of X-particles i.e. 
if feedback is present, the properties of nz can be derived formally as a function of n and 
then the properties of n can be derived self-consistently trough ( 1721) . as done in subsection 

(HE]). 

IV. COMMENTS AND CONCLUSIONS 

In this paper we have analyzed general stochastic birth and death models that include 
delay. We have presented three different methods that together constitute a general toolbox 
to study stochastic models including delay. 
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Poissonian 
Hopf 



Super-Poissonian 



Sub-Poissonian^ 



FIG. 3: Relative size of the variance respect mean value for the number of particles, in the r — e 
plane, for creation with constant delay and a negative feedback given by a creation rate of the form 
i+(e<j))' 2 ( no * e th & t e is the strength of the negative feedback). The "Poissonian line", a 2 [n\ = (n), 
obtained through the approximation (|68p . marks the transition from sub-Poissonian to Super- 
Poissonian fluctuations, while the Hopf line marks the Hopft transition into oscillatory behavior in 
the deterministic system. Parameters values are: c = 1, D = 7 = 1. (Online version in colour.) 



In sub-section (III Al) we have shown that when the creation rate is independent of the 
state of the system (no feedback) and the initiation of the delayed degradation and the 
instantaneous degradation are first order reactions (rate not depending on the state of the 
system), the process can be solved fully in an exact fashion for general distributions of 
delay, showing always Poissonian character and a monotonically decreasing time correlation 
function given by (JS}. 

In sub-sections ( III Bj) . ( Ill Cj) we have considered a more general process with delay in 
the degradation step, allowing the initiation of the delay degradation and the instantaneous 
degradation to be higher order reactions, as well as the presence of feedback in the creation 
rate. The method allows to reduce the system to a Markovian one, where usual techniques 
can be used. Explicit expressions for the time correlation for general delay distributions 
were obtained. In this case the correlation might be non-monotonic, if feedback is present, 
but typically decreases monotonically. 

Section f illip shows that when the delay appears in the creation reaction and feedback 
is present, the delay typically has more dramatic consequences. In the case of fixed delay, 
it is shown that for negative feedback, the fluctuations are amplified as the delay increases, 
going beyond the level found when no feedback is present, and the time correlation function 
becomes oscillatory, alternating positive and negative values at approximately multiples of 
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the delay. In the positive feedback case, again for fixed delay, the fluctuations are reduced 
with increased delay and the time correlation function remains always positive. 



V. APPENDIX: CALCULATION OF P(n,t) IN THE SIMPLE CASE OF DE- 
LAYED DEGRADATION 

We start by considering the case n — 0. For the sake of simplicity, we focus on the case 
with creation rate, C, independent of time, but the generalization to time-dependent C is 
straightforward. Since the time origin is taken at t — 0, the probability of observing zero 
particles at time t > is equal to the following limit: 

M-1 

P(0, t) = lim TT [1 — CAt + CAtF{t - U) + o(At)} , (75) 

i=0 

with At = j- playing the role of a small time- increment and ti = iAt. This expression 
follows from the fact that, in order to find the system with zero particles at time t, in every 
previous infinitesimal time interval (£' G [ti, t i+1 ),i — 0, . . . , M — 1) one of the following two 
(incompatible) events must take place: either a particle is not created (probability 1 — CAt) 
or a particle is created with a lifetime smaller that t — U (probability CAtF(t — ti)). We 
now have: 

M-1 r t 

logP(0,t)= lim V \-CF(t-ti) + o(At) At = -C / dt'F(t-t'), (76) 

M->oo ^— ' L J In 

with F(t) = 1 - F(t), so we find 

P(0,f)=e-° ■£<*'*■(*-*'). (77) 
Following a similar line of reasoning, P(n, t) can be computed as: 

M-1 M-1 M-1 n 

P(n,t) = Jirn^ E Ilt '^ (*"**«)] II [l-CAtP(t-tg 

ii=0 i 2 =H+l i n =j„_l+l J=l 0<3<M-1 

(78) 

This expression results from the consideration of choosing the times (t^, . . . ,t in ) at which 
the n particles are created and survive up to t. The Z-th particle is created with probability 
CAt and survives up to t with probability F (t — t^). The other factor comes from the fact 
that at the other time intervals either a particle is not created or it is created but dies before 
t. 
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Using 



lim 

M-»oo 



] ] [l - CAtF(t - U 



0<j<M-l 



-Cf dt'F{t-t>) 



(79) 



and replacing the sums by integrals in the limit M — > oo 

dtxCF(t - h) / dt 2 CF{t - t 2 ) . . . / dt n CF{t -t n ) = — 



tn-l 



nl 



dt'F(t - 1) 



we finally obtain: 



p^ t ) =e -Cj^t'F(t-t')_ 



Cdt'F{t-t') 



10) 



(81) 



that is, a Poisson distribution with average (n(t)) = C j^dt'F{t — t'). In the steady state 
(found as the limit t — > oo), the average becomes {n(t)) = C(t). Remarkably, this Poissonian 
character is completely independent of the form of the delay distribution. As commented 
above, this result can be easily generalized to the case in which the creation rate depends 
on time, C — )■ C(t), obtaining again a Poisson distribution with average J* dt'C{t')F{t — t'). 



VI. APPENDIX: DERIVATION OF THE MASTER EQUATION IN A SYSTEM 
WITH DELAY 

Here we derive the master equation of the process (153|) . We consider first the case of fixed 
delay r. We start with the following identity: 

P(n,t+A) = ^2P(n,t+A;n',t) = P(n,t+A;n+l,t)+P(n,t+A;n-l,t)+P(n,t+A;n,t)+o(A). 

n' 

(82) 

It is immediate to see that P(n, t + A; n+ 1, t) = ^{n + l)AP(n + 1, t). In the case of fixed 
delay, the second sum can be evaluated introducing a tree-times probability as: 

P(n,t + A;n - l,t) = P(n, t + A; n - 1, t; n', t - r) 

n' 

= p ( n i t + A\n- 1, t; n', t - r)P(n', t-r;n-l,t). (83) 

Now, P(n, t + A|n — 1, t; n', t — r) = C(n')A + o(A). Expanding in a similar way the term 
P(n,t + A;n,t), and taking the limit A — > 0, we can obtain the master equation of the 
process: 



= (E-l^nPin^ + iE- 1 -!) 



^TC(n')P(n',t-T;n,t) 



n'=0 
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In the case of distributed delay, we start considering a discrete distribution of delays i.e. 
t = n, . . . ,Tm with corresponding probabilities /(ti), . . . /(tm). The continuum limit can 
then be obtained making M — > oo. The creation term in (j82p can be written as: 

P(n,t + A;n — l,t) — P(n,t + A;n — l,t;ni,t — Ti; . . . ;n M ,t — r M ) (85) 

= 2J P(n,t + A|n- l,t;ni,t - . . r,n M ,T m )P{ni,t - ti; . . . ;n M ,t - T M ;n - l,t). 

n-i,...,n M 

Now, P(n, t + A|n — 1, t; ni, t — r 1; n M , r m ) = 1 C(rij)/(rj)A + o(A), that is, the 
probability that a particle started its creation at time t — Ti with a creation time equal to 
Tj. Replacing in the previous equation and performing the appropriate sums we obtain: 

M 

P{n, t + A;n-l,t) = ^2Y^ C ( n ')f( T i)P(n', t - r<; n - 1, t) A + o(A) (86) 

n' i=l 

that in the continuum limit reduces to J2 n ' Io° d T C(n')f(T)P(n',t — T;n — l,t). Considering 
in a similar way the other terms in fl82|) and taking the limit A — > one can obtain the 
master equation for distributed delay (!56|) . 
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